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Abstract. In this note we show that the non-symmetric version of the classical Nitsche's method 
for the weak imposition of boundary conditions is stable without penalty term. We prove optimal 
i? 1 -error estimates and L 2 -error estimates that are suboptimal with half an order in h. Both the 
pure diffusion and the convection-diffusion problems are discussed. 

1. Introduction. In his seminal paper from 1971, [TS], Nitsche proposed a con- 
sistent penalty method for the weak imposition of boundary conditions. The formula- 
tion proposed was symmetric so as to reflect the symmetry of the underlying Poisson 
problem. Stability was obtained thanks to a penalty term, with a penalty parameter 
that must satisfy a lower bound to ensure coercivity. 

A non-symmetric version of Nitsche's method was later proposed by Freund and 
Stenberg [9] and it was noted that this method did not need the lower bound for 
stability. The penalty term however could not be omitted, since coercivity fails, and 
error estimates degenerate as the penalty parameter goes to zero. The non-symmetric 
version of Nitsche's method was then proposed as a discontinuous Galerkin (DG) 
method by Oden et ai., |16] and it was proven that the non-symmetric version was 
stable for polynomial orders k > 2, by Girault et al. [17] and Larson and Niklasson 
[2]. In [H] stability for the penalty free case is proved using an inf-sup argument 
that relies on the important number of degree's of freedom available in high order 
DG-methods. 

To the best of our knowledge no similar results have been proven for the non- 
symmetric version of Nitsche's method for the imposition of boundary conditions when 
continuous approximation spaces are used. Indeed in this case the DG-analysis does 
not work since polynomials may not be chosen independently on different elements 
because of the continuity constraints. Weak impositition of boundary conditions has 
been advocated by Hughes et al. for turbulence computations of LES-type in pQ. 
They showed that the mean flow in the boundary layer was more accurately captured 
using weakly rather than strongly imposed boundary conditions. They also noted 
that the non-symmetric Nitsche's method appears stable without penalty |13| . 

In applications there is interest in reducing the number of free parameters used 
without increasing the number of degrees of freedom needed for the coupling, see 
[TU] for a discussion. From this point of view a penalty free Nitsche method is a 
welcome addition to the computational toolbox, in particular for flow problems where 
the system matrix is non-symmetric anyway, because of the convection terms. It has 
no penalty parameter and does not make use of Lagrange multipliers. 

Numerical evidence also suggests that the unpenalized non-symmetric Nitsche 
type method has some further interesting properties. When using iterative solu- 
tion methods in domain decomposition it has been shown to have more favorable 
convergence properties compared to the symmetric method [5]. For the solution of 
Cauchy-type inverse problem using steepest descent type algorithms it has been shown 
numerically to have superior convergence properties in the initial phase of the itera- 
tions compared to the symmetric version or strongly imposed conditions, in spite of 
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the lack of dual consistency. 

In view of this the question naturally arises if the penalty free method is sound, 
or if it can fail under unfortunate circumstances. 

In this paper we prove for the Poisson problem that the non-symmetric Nitsche's 
method is indeed stable and optimally convergent in the iJ 1 -norm for polynomial 
orders k > 1. We also show that in this case, the convergence rate of the error in 
the L 2 -norm is suboptimal with only half a power of h. Hence the non-optimality 
for the non-symmetric Nitsche's method for continuous Galerkin methods is not as 
important as for DG-mcthods (see [TB] and [TT] for numerical evidence of the sub- 
optimal behavior in this latter case). 

We then show how the results may be applied in the case of convection-diffusion 
equations, considering first the Streamline-diffusion method and then outlining how 
the results may be extended to the case of the Continuous interior penalty method. 

Nitsche's method however has some stabilizing properties of its own, in particular 
for outflow layers, this phenomenon was analyzed in 1SJ and illustrated herein with 
a numerical example. This makes the non-symmetric Nitsche's method an appealing, 
parameter free, method for flow problems where the system matrix is non-symmetric 
and the use of stabilized methods usually also results in the loss of half a power of 
h. It should be noted however that the smallest error in the L 2 -norm is obtained 
with the formulation using penalty on the boundary, as illustrated in the numerical 
section. So we do not claim that the penalty free method is the most accurate. 

We only prove the result in the case of the imposition of boundary conditions but 
the extensions of the results to the domain decomposition case of [5] or the fictitious 
domain method of [1] are straightforward using similar techniques as below. Also 
note that since the main aim of the present paper is the study of weak imposition of 
boundary conditions, we will assume that the reader has basic understanding of the 
techniques for analyzing stabilized finite element methods and some arguments are 
only sketched. 

For the sake of clarity, we first prove the main result on the pure diffusion prob- 
lem and then discuss the extension of our result to the case of convection- diffusion 
problems. The paper is ended with some numerical examples. 

2. The pure diffusion problem. Let ft be a bounded domain in M 2 , with 
polygonal boundary <9f2. Wherever iJ 2 -regularity of the exact solution is needed we 
also assume that fl is convex. Let {r^ji denote the faces of the polygonal such that 
<90 = U^. The Poisson equation that we propose as a model problem is given by 



—Au — / in fi, 
u = g on dCl, 



(2.1) 



where / G L 2 {Q) and g G iJ 1/2 (<9ft) or g G H^ 2 (dQ). 

We have the following weak formulation: find u G V g such that 



a(u,v) = (/,«)n, Vv eVq, 



(2.2) 



where (x, y)o denotes the L 2 -scalar product over O, 

V g :={veH 1 (Q):v\ dn = g} 



and 



a(u,v) := (Vu, Vv)q. 
2 



This problem is well-posed by the Lax-Milgram's lemma, using the standard ar- 
guments to account for non-homogeneous boundary conditions. The i/^-stability, 
IMIff^fi) < Ci?i(||/|| + llfflliri/s^fi)) holds and under the assumptions on ft, f and g 
there holds \\u\\ H 2 (n) < C m {\\f\\ + ||ff|| ff 3/2 (an) ). Here we let ||x|| := \\x\\ L 2 {n) . Below 
C will be used as a generic constant that may change at each occasion, is independent 
on h, but not necessarily of the local mesh geometry. We will also use the notation 
a < b for a < Cb. 

3. The finite element formulation. Let {Th} denote a family of quasi uniform 
and shape regular triangulations fitted to f2, indexed by the mesh-parameter h. (It 
is straightforward to lift the quasi uniformity assumption, at the expense of some 
standard technicalities and readability.) The triangles of Th will be denoted K and 

o 

their diameter Hk '■= diam(AT). The interior of a set P will be denoted P. For a 
given Th the mesh-parameter is determined by h := vasxKeT h hn- Shape regularity is 
expressed by the existence of a constant c p G K for the family of triangulations such 
that, with pk the radius of the largest ball inscribed in an element K, there holds, 

~< Cp ,VKeT h . 

PK 

For technical reasons, and to avoid the treatment of special cases, we assume that for 
all i, Ti contains no less than five element faces. 

We introduce the standard finite element space of continuous piece wise polyno- 
mial functions 

V* := {v h € H l (fl) : v h \ K E P k (K), VK E T h }, k > 1, 

where ¥k(K) denotes the space of polynomials of degree less than or equal to k on 
the clement K. The finite clement formulation that we consider then takes the form, 
find Uh € such that 

a-h(uh, Vh) = (/, Vh)n + (g, V«ft • n) gn \fv h E V£, (3.1) 

where (x, y)g a denotes the L 2 -scalar product over the boundary of fl and 

a h (u h ,v h ) := a(u h ,v h ) - (Vuh ■ n,v h ) 9n + (u h ,Vv h ■ n) 9n . (3.2) 

Note that in the classical non-symmetric Nitsche's method we also add a penalty term 
of the form 

^ {ih-KUh, v h) dnn9K (3.3) 
K 

and modify the second term of the right hand side accordingly 

Y.i^^Vh + Vv h ■ n) dnndK . 

K 

The key observation of the present work is that the penalty parameter 7 may be 
chosen to be zero without loss of neither stability nor accuracy. 



Inserting the exact solution u into the foundation (3.1) and integrating by parts 
immediately leads to the following consistency relation. 

3 



Lemma 3.1. If u is the solution of (2.1) and uh is the solution of (3.1) then 
there holds 

a h (u - u h ,v h ) = 0. 

For future reference we here recall the classical trace and inverse inequalities satisfied 
by the spaces V*. 

Lemma 3.2. (Trace inequality) There exists Ct € M such that for all Vh € Pfe(if) 
and for all K G 7ft there holds 

_i i 
\\vh\\L 2 (dK) < C T (h K 2 \\v h \\ L 2( K) + h 2 K \\^v h \\ L 2 {K) ). 



Lemma 3.3. (Inverse inequality) There exists Cj € K such that for allvh € ¥k{K) 
and for all K £ % there holds 



\\^Vh\\L^(K) < Clhjl\\Vh\\L 2 {K)- 



4. Stability. The non-symmetric Nitsche's method is positive and testing with 
Vh = Uh immediately gives control of the i/^-seminorm of Uh- In order for the formu- 
lation to be well-posed this is not sufficient. Indeed well-posedness is a consequence 
of the Poincare inequality that holds provided we have sufficient control of the trace 



of Uh on dfl. This is the role of the penalty term (3.3), it ensures that the following 
Poincare inequality is satisfied 

IKH < C P \\u h \\ llh , where \\u h \\l h := \\Vu h \\ 2 + |K|||, Mn 

with 

IKHi,Mfi : = J2 ( h ~K u h,u h ) dnn9K . 

K 

Since we have omitted the penalty term, boundary control of Uh is not an immediate 
consequence of testing with Vh = u/i- What we will show below is that control of 
the boundary term can be recovered by proving an inf-sup condition. Indeed the 
non-symmetric Nitsche's method can be interpreted as a Lagrange multiplier method 
where the Lagrange multiplier A?, has been replaced by the normal gradient of the 
solution: Vu/j • n. This interpretation of the Nitsche's method was originally proposed 
in |20j . however without considering the inf-sup condition. In the DG-framework 
it was considered in [7], where equivalence was shown between a certain Lagrange- 
multiplier method and a certain DG-method. When Lagrange-multipliers are used to 
impose continuity, the system has a saddle point structure and the inf-sup condition is 
the standard way of proving well-posedness. Here we will follow a similar procedure, 
the only difference is that the solution space and the multiplier space are strongly 
coupled, since the latter consists simply of the normal gradients of the former. A 
key result is given in the following lemma where we construct a function in the test 
space that will allow us to control certain averages of the solution on the boundary. 
To this end regroup the boundary elements, i.e. the elements with either a face or a 
vertex on the boundary, in (closed) patches Pj, with boundary dPj, j — l...Np. Let 
Fj := dPj n <9f2. We assume that the Pj are designed such that each Fj has at least 



four inner nodes (this is strictly necessary only if both end vertices of Pj belong to 
corner elements with all their vertices on the boundary). Under our assumptions on 
the mesh, every contains at least one patch Pj and there exists ci , C2 such that for 
all j 

c\h < meas(Fj) < c%h. (4.1) 

The average value of a function v over Fj will be denoted by vP . 

Lemma 4.1. For any given vector (rj) € M. Np there exists ip r € such that 
for all 1 < j < Np there holds 

meas(Fj)^ 1 / \7tp r ■ n ds = rj (4-2) 

J Ft 



\ 1/2 

Proof. We first construct a function tpj taking the value 1 in the interior nodes 
of dfl n dPj and zero elsewhere. Fix j and let (pj £ be defined, in each vertex 
Xi € %, by 

{0 for Xi G K such that K has three vertices on dil; 
for Xi € O \ ; 
1 fOr lEj G Pj. 

Let 



meas(Pj) 1 / V<^j • n ds 
Jf, 

and define the normalised function tpj by 



<Pj : = S j Vj- 



This quantity is well defined thanks to the following lower bound that holds uniformly 
in j and h 

C* H < Ejh. 

The constant C~, only depends on the local geometry of the patches Pj . By definition 
there holds 

meas(Fj) -1 / Vifj ■ n ds = 1 (4.4) 

J Fi 



and using the standard inverse inequality (Lemma 3.3) 



IV^II < dh-^W^U^ < C / /i- 1 5- 1 meas(P i ) 1/2 < CjC^h. (4.5) 
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Now defining 



N P 

3=i 



we immediately see that condition (4.2) is satisfied by equation (4.4). The upper 
bound (4.3) follows from (4.5), the relation (4.1) and using that 

Np 

■= \\h~ir 112 



l^lli h,dn 



N P 



N P 



With the help of this technical lemma it is straightforward to prove the inf-sup 
condition for the formulation (3.1). 

Theorem 4.2. There exists c s > such that for all functions v^, £ Vt there 
holds 



C s \\Vh\\l,h < SUp 



a h (v h ,w h ) 



Uh£V£ \\ W h\\l,h 

Proof. Recall that 

a h (v h ,w h ) = (Vv h , Vw k )si - (Vv h ■ n,w h ) dn + (v hl Vw h ■ n) m 
Taking w h = v h gives 

a h (v h ,v h ) = \\VvhW 2 . 
To recover control over the boundary integral we let 

rj = h~ l vi := /i _1 meas(F :? ) _1 / Vh ds 

Jf 3 

in the construction of ip r in Lemma |4.1| and note that 

N P 

>V<Pr -n) en = J2 (ll^" 1/2 ^ll^(F 3 .) + (K ~ V<p r ■ n) F , 



(Vh, 



Using standard approximation, 

IK - v 3 \\ L 2 {F . } < h\\Wv h x n\\ L 2 {Fj) , 



(4.6) 



(4.7) 



and by the trace and inverse inequalities of Lemma 3.2 and Lemma 3.3 we have 

((vh - v r ),W<f r ■ n) p < C T (1 + C/)||Vw, l || L 2(p j) ||V(^ r || L 2 ( p j) . 
Moreover since by Cauchy-Schwarz inequality and the trace inequality, 
\(Vv h , Vwh)n - (Vv h ■ n,w h ) gn \ < \\Vvh\\\\w h \\i <h 
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we deduce using the stability (4.3) that 



N P 



a h (v h , ip r ) > E \\h- 1/2 v J \\h {Fj) - C\\Vv h \\\\<pr\\i, h 



3=1 



Np 



> E ll^^lli^F,) - aiiv^n E II^^H^CF,) 



3=1 



We now fix = Vf l + nip r and note that 



E 



1/2 



nJ'll?. 



o h (« h ,«; h ) > ||V^|| 2 + r ? EH /l " 1/V l!i 2 (^) 



E 

3=1 



E 

v3'=l 



1/2 



anv^iir? Eii^ 1/V ii^w- 



rj II ?. 



Np 



> (1 - £ )||V^|| 2 + ^(1 - <7 s V(4e)) E ll^^'lli^,). ( 4 - 8 ) 

3=1 

It follows, using once again the approximation properties of the L 2 -projection on the 



piece wise constants (4.7), that for any e < 1 we may take rj sufficiently small so that 



there exists <v e such that 



Np 



Cr,,e\\v h \\l h < CCr,, e \ \\ Vu h \\ 2 + E II h 1/2 V 3 Wh^ ] < a h (v h ,W h ). 



We may conclude by noting that by (4.3), our choice of rj and the stability of the 
L 2 -projection on piece wise constants there holds 



\wh\\l,h < \\vh\\l,h + Tl\\<Pr\\l,h < Cjj||«fc||l,fc- 



(4.9) 



□ 



5. A priori error estimates. The stability estimate proved in the previous 
section together with the Galerkin orthogonality of Lemma [3~T] leads to error estimates 
in the || • || i ^ norm in a straightforward manner. First we will prove an auxiliary lemma 
for the continuity of a/j(-, •). To this end we introduce the norm 

INI* := IMIi.fr + \\h%Vu • n\\ L 2 {dn y 



Lemma 5.1. Let u e H 2 (fl) + V£ and Vh € V*. Then the bilinear form a h (-, •) 
defined by (3.2) satisfies 



ah{u,v h ) < C\\u\\*\\v h \\ l h . 



Proof. The result is immediate by application of the Cauchy-Schwarz inequality 
and the trace inequality of Lemma |3.2| □ 
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PROPOSITION 5.2. Let u e H k+1 (fl) be the solution of ( |2.1[ ) and u h the solution 
of (3.1). Then there holds 



\\u - u h \\i,h < Ch k \u\ H k+i {n) . 

Proof. Let i^ z u denote the Scott-Zhang interpolant of u jTHj. Using the approxi- 
mation properties of the interpolant it is straightforward to show that 



\u - is Z u\\iJi + \\u - is Z u\\ 



< h k \u\ H k+i {n) . 



We therefore use the triangle inequality to obtain 

||« - u h \\i,h < \\u - isz u h,h + \\u h - isz u \\i-h> 

where only the second term needs to be bounded. To this end we apply the result of 
Theorem |4.2| followed by the consistency of Lemma |3.1| 



\u h - Hz u \\i,h < sup 



ah{uh - it z u,w h ) 



w h £V« 



\ w h\\l,h 



sup 



a h (u - i$ z u,Wh) 



\\ w h\\l,h 

By the continuity of Lemma 5.1 and the approximation properties of i^ z u we conclude 



C s \\Uh - hz u \\i,h 



< \\u-i k u\\ 



□ 



For DG-mcthods it is well-known that the non-symmetric version may suffer from 
suboptimality in the convergence of the error in the L 2 -norm due to the lack of adjoint 
consistency. This is true also for the non-symmetric Nitsche's method, however since 
the method is used on the scale of the domain and not of the element the suboptimality 
may be reduced to hfl as we prove below. 



Proposition 5.3. Let u € H k+1 (Q) be the solution of (|2.1[) and u h the solution 



of (3.1). Then 



\\u-u h \\ < Ch k+1 i\u\ Hk +i {n) 



Proof. Let z satisfy the adjoint problem 

J — Az = u — Uh in f2, 
1 z = on dfl. 

Under the assumptions on f2 we know that ||z||ff2(n) < Cr2\\u — Uh\\- It follows that 

|| ti - u h \\ 2 = (u - u h , -Az)q = (V(u - u h ),Vz) n - (u- u h ,Vz ■ n) dn 

= a h (u - u h , z) + 2(u- u h ,Vz ■ n) m . 

By Lemma |3.1| and a continuity argument similar to that of Lemma |5.1| using that 
(z — i\ z z)\q£i = 0, it follows that 

a h (u - u h , z) = a h (u - u h ,z - i\ z z) 

= (V(w - u h ),W(z - il z z)) n - (u - u h , V(z - i\ z z) ■ n) gn 
< \\u - u h \\i,h\\z ~ 4z 2 ll* 

< h\\u - Uh\\\,h\z\H 2 (n)- (5-1) 
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We also have, using the following global trace inequality 



||Vz ■ n\\ L 2 (m) < ||z|| H 2 (0) , 

that 

| (u - u h ,Vz ■ n) an | < h 1/2 \\u - u h \\i hm \\z\\ H 2 {n y (5.2) 



Collecting the inequalities (5.1) and ( |5.2[ ) we arrive at the estimate 
\\u-u h \\ 2 < (h + h 1/2 )h k \u\ H k+i m \\z\\ H 2 m 

and we conclude by applying the regularity estimate ||z||if2(n) < Cr2\\u — Uh\\- □ 

6. The convection— diffusion problem. Since the method we discuss leads to 
a non-symmetric system matrix the main interest of the method is for solving flow 
problems where an advection term makes the problem non-symmetric anyway. Note 
that there appears to be no analysis that is robust with respect to the Peclet number, 
even in the case of the non-symmetric discontinuous Galerkin method. 

We will therefore now show how the above analysis can be extended to the 
case of convection-diffusion equations yielding optimal stability and accuracy both 
in the convection and the diffusion dominated regime. We will consider the following 
convection-diffusion-reaction equation: 

(ju + fi- Vu-eAu = f in O, (6.1) 

and homogeneous Dirichlet boundary conditions. We assume that f3 £ [W^ (^)] 2 , 
cr G R, 

o - * V • p > c CT > 

and e € R + . In this case the formulation writes: find E 14 such that 

Ah(uh, v h ) := (au h + /3 ■ Vu h , v h )u - (P ■ n, u hl v h ) m - 

+ ea h (u h ,v h ) = (f,Vh)n, Vti/, e 14 (6.2) 

where := {x g d£l : ±/3 • n > 0}. First note that the positivity of the form now 
writes 

A h {u h ,u h ) > ^|!|/?-n|5 U/l ||? )0 + ||e5Vu /l || 2 , (6.3) 

hence provided \(3-n\ > on some portion of the boundary with non-zero measure the 
matrix is invertible. In the following we assume that this is the case, but we do not 
assume that \j3-n\ > everywhere on dtt. To prove optimal error estimates in general 
we require stronger stability results of the type proved above to hold. It appears 
difficult to prove these stronger results independently of the flow regime. Indeed it is 
convenient to characterize the flow using the local Peclet number: 

Pe := m. 



If Pe < 1 the flow is said to be diffusion dominated and if Pe > 1 we say that it is 
convection dominated. We will now treat these two cases separately. 
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In view of the equality (6.3) we introduce the following strengthened norm 



IKIIiA/3 '■= e hh\\l h + \\\\P ■ n\^v h \\ 2 m . 

This norm is suitable in the diffusion dominated regime, but will be modified by the 
introduction of stabilization when the convection dominated regime is considered. 

6.1. Diffusion dominated regime Pe < 1. In this case we may prove an inf- 
sup condition similar to that of Theorem |4.2| For simplicity we assume that a = 0. 



Proposition 6.1. (Inf-sup for convection-diffusion, Pe < 1.) For all functions 
Vh G Vfi there holds 

A h (v h ,w h ) 

Cs\\vh\\i,h,p < sup -J, n . (6.4) 

Wh ev£ \\ w h\\i,h,p 



Clearly, compared to the proof of Theorem 4.2 we only need to show how to handle 
the term 

(f3 ■ Vv h , (fr)n - (P ■ nv h , tp r ) Bn - ■ 
The necessary bound on this term is given in the following Lemma. 



Lemma 6.2. Let ip r be the function of Lemma J^.l with r chosen as in ( 4.6 1 . 
Then for Pe < 1 there holds for all \i > 

(f3 ■ Vv h , rnp r )u - (/3-nVh, rjip r ) gn - 

< ^e\\Vv h \\ 2 + HI/? ■ n\^v h \\ 2 LHm) ) + C 2 d (^)-We\\v h \\\ h dn . 

Proof. Let 

(P ■ Vv h ,r)<p r )u - (0 • nv h ,nLp r ) gn _ = T x + T 2 . 
By the definition of the Peclet number and the Cauchy-Schwarz inequality, we have 

Tx < Peel\\Vv h \\riet\\h- l (pr\\. 



From the construction of ip r , a scaling argument, the stability (4.3) and the choice of 



r (4.6) we deduce that 

II^Vll<l|v^<c s K||i A9n . 

Using the arithmetic-geometric inequality we have 

Tx < /lellV^H 2 + C 2 d {^)- l Pe 2 ife\\v h \\\ hda . 
For Ti we have using a Cauchy-Schwarz inequality, the definition of the Peclet number 



and the stability (4.3) 



T 2 < \\\(3 ■ n\2v h \\ L 2 {dn) Pe2ne^\\ip r \\^ htdn < C d \\\0 ■n\^v h \\ L 2 {dn) T 1 e^\\v h \\i hdn . 
We apply the arithmetic-geometric inequality once again to conclude. □ 
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Proof. (Proposition 6.1) The inf-sup stability (6.4) now follows by taking Wh '■— 
vi l +r)ip r and proceeding as in equation (4.8 1 using (6.3 ) and Lemma 6.2 in the following 
fashion 



A h {v h , v h + wv) > (1 - e - p)e\\ Vv h \\ 2 + {\~ A»)|| \P • n\^v h \\l 2(n) 

+ 77(1 - C^/(4e) - C^/(4M))£||^|||, ft ,aa- 

We may now choose e = 1/4 and /Lt = 1/4 and then small enough so that positivity 
is ensured. Then 

A h (v hl v h + rjipr) > C v \\vh\\i tht p- 



We conclude as in Theorem 4.2 but now using the norm || • 



\\w h \\lM,p < \\Vh\\l,h,P + 1l\\<fr\\l,h,f3 < \\vh\\l,h,p + vC\\vh\\i,h,P + vlWP ■ ra| 2 Vr||i 2 (an) 

□ 

Proceeding as in Proposition |5.2[ this leads to optimal a priori estimates in the 
norm II • \U h for Pe < 1. 



Proposition 6.3. Let u € i? fc+1 (fi) &e t/ie solution of (|6.1[) and it/, ifte solution 



of (6.2) and assume that Pe < 1. TTien 



Proof. As in the proof of Proposition |5 . 2| we arrive at the following representation 
of the discrete error 

n - k n ^ A h {u h - is Z u,w h ) A h (u-i% z u,w h ) 
c s \\uh — ia Z u\\i,h,p < sup n n = sup 



tu h evj \\ w h\\i,h,p w k ev* \\Wh\\l,h,P 



By the continuity of Lemma 5.1 and an integration by parts in the convective term 
we obtain 

A h (u h - is Z u,w h ) < e\\u - is Z u\\*\\ w h\\i,h 

+ (u- is Z u, P ■ Vw h ) Q + (f3 ■ n(u - i% z u), w h ) gQ+ ) 
^ e(||tt-isHI* + Pe\\h~ l {u - ^ Z u)|| + Pe\\u - it z u\\x <Km )\\w h \\i,h,p- 

As a consequence 

£\Wh - itz u \\iJi < \\uh - isz u \\i,h,P 

< c- l e{\\u - i* z u||, + Pe\\h-\u - i k sz u)\\ + Pe\\u - &u|| i iMtJ )- 

The claim follows by dividing through by e, using approximation and the assumption 
Pe < 1. □ 
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6.2. Convection dominated regime: the Streamline diffusion mehod. 

In the convection dominated regime, when Pe > 1, we need to add some stabilization 
in order to obtain a robust scheme. We will here first consider the simple case of 
Streamline-diffusion (SD) stabilization and assuming a = 0. In the next section the 
results will be extended to include the Continuous interior penalty (CIP) method. 
The formulation now takes the form: find Uh € Vt such that 

A S d{uh, v h ) := (P ■ Vlffc, v h + SP- Wv h ) n -^(eAu^SP ■ Vv h ) K - (P ■ n u h ,v h ) du - 

K 

+ sa h {u h ,v h ) = {f,v h + SP-Vv h ) n , V» 4 G^ (6.5) 

where S = r )soh/\P\ when Pe > 1 and 5 = otherwise. At high Peclet numbers, the 
enhanced robustness of the stabilized method allows us to work in the stronger norm 
IH^HIm defined by 

IKIHL := • Vu h f + h\p ■ n\^u h \\l, {m) +e\\Wu h f. (6.6) 



We will also use the weaker form |||u^|)| 2 defined by (6.6) with (5 = and for the 
convergence analysis we introduce the norm 

K 



Testing the formulation (6.5) with = Uh yields the positivity 

c|||«h|||ft i( 5 < A SD (uh,u h ) (6.7) 

in the standard way using an element wise inverse inequality to absorb the second 
order term, i.e. 

J2(eAu h ,6p ■ Vu h ) K < \chsDPe- 1/2 \\e^Vu h \\ 2 + h\6^p ■ Wu h \\ 2 . 

K 

Clearly for jsd < 1/(C|) stability holds for Pe > 1. 

Unfortunately the norms proposed above seem too weak to allow for optimal error 
estimates. Indeed, since we do not control all of ||ufc||i,/„ for general u € H 2 + V h l , 
Vh € there does not hold Agjj(u,Vh) < IIMIMIkfilllM' ( c -f- Lemma 5.1) unless an 
assumption on the boundary velocity such as \P ■ n\h > e is made. It also appears to 
be difficult to obtain an inf-sup condition similar to (6.4) in the high Peclet regime. 

We therefore use another technique to prove optimal convergence directly. The 
idea is to construct an interpolation operator ttqu, such that the interpolation error 
u — ngu satisfies the continuity estimate: 

A S d{u - 7T9U,v h ) < \\\u - 7r a u|||*||Mlk<5- (6.8) 

Assume that we have an interpolation operator ng : i/ 1 (SI) such that the 

following hypothesis are satisfied. 
(HI) Approximation, 

\\TT d u - u\\ + h\\\7(ir 9 u -u)\\< Ch k+1 \u\ Hk+ i {n) . (6.9) 
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(H2) Normal gradient, 



/ \7(ttqu — u) ■ n ds = 0, i = 1 . . . Np, 



(6.10) 



where Fi are the boundary segments introduced in Section [4j 
Under assumptions (HI) and (H2), we may prove the optimal convergence of the 
SD-method. 



Proposition 6.4. Let u € H k+1 (Q) be the solution of (|6.1[) and Uh the solution 



of (6.5). Assume that there exists ttqu € V ft satisfying (H.1) and (H.2). Then 
|||«-«h||lM ~ /i fe+ ^(l + ^e~5>)| u | fffc+1(n) . 



Proof. It follows from the approximation properties of 7ra that 
III" - * a u|||, < \\f3\\ih k+ Hl + Pe-^ 2 )\u\ H H +Ha) . 



We now need to prove the continuity (6.8). Note that 



Asd{u - n d u,v h ) = V(u - ir d u) + S *(u- ir d u),Sif3 ■ ^Jv h ) K 

- ^2(S^eA(u - n d u),5^(3 ■ Vv h ) K + ([3 -n(u- Tr 9 u),v h ) gn+ + ea h (u - n 9 u, v h ) 



K 



< 



\\\u - 7rau|||* ||Mlk<5 + ea h (u - TT d U,V h ) . 



Consider now the term I\. We will prove the continuity 

ea h (u-Tr d u,v h ) < e*\\u-ndu\\*ez\\vh\\i,h < \\\u ~ ^du\\\*\\\ v h\\\h,s (6-11) 

Using Cauchy-Schwarz inequality and a trace inequality we show the continuity of the 
first and last term. 



■ (V(u - 7r 9 u), V^) n - e (V(u - tt 9 u) ■ n, v h ) m + e (Vv h -n,(u- ir d u)) 



'an 



< || u - 7rau||*|||vft||U,o - £ (V(u - ttou) ■ n,v h ) 



For the remaining term we must exploit the orthogonality property ( 6.10[ ) of ngu on 
the boundary. Indeed by decomposing the boundary integral on the Np subdomains 
Fi we have, denoting by v l h the average of Vh over the boundary segment Fi- 



N P 



e (V(u - irgu) ■ n,v h ) gn = (V(it - ngu) -n,v h - v l h ) Fi 



N P 



< 



£^2\\^( u - nan) ■ n\\ L 2 {Fz) \\v h - v l h \\ L 2 (Fz) 



< e 5 1 1 V (u - ir d u) ■ n 1 1 _ i X8n e 5 1 1 Vv h | | 

< e^||«-7rsu||«e2||u h || 1)fe . 

13 



Where we used the approximation properties of the local average and a trace inequal- 
ity. Collecting the above estimates and noting that 

ei\\u - n d u\\* < \\\u - 7r a u|||», e^\\v h \\i. h < ||K||U,o 



concludes the proof of (6.8). 



Using the positivity (6.7), and the consistency of the method wc have, setting 



e;,. := Uh — ttqu, and using that Pe > 1 

IIKIII/U = A S D(eh,e h ) = A S d{u - ndU,e h ) < \\\u - 7rait|||*|||e h ||| fc)5 

<h k+ i \\0\\L(1 + Pe-*)|«|H* + i (n) ||| e/l ||| M . 

□ 

We end this section by the following Lemma establishing the existence of the inter- 
polation ng with the required properties. 

Lemma 6.5. The interpolation operator ttq : H 1 ^) i— >• satisfying the proper- 
ties (H.1 ) and (H.2 ) exists. 

Proof. Let ttqu := i^ z u + ip r where ip r is the function of Lemma 4.1 with the rj 
chosen such that 

Tj = Vu • n — Vig Z u • n . 
Clearly by construction there holds 

/ (W*. • n - V. ■ n) ds = / (V&« • n + V„ r ■ n - V. ■ n) d, 

JFi JFi 

= / (Vig Z w • n + Ti - Vu ■ n) ds = 0. 

To prove the approximation results we decompose the error 

\\u - n d u\\ < \\u - is Z u\\ + \\is Z u - n d u\\ < Ch k+1 \u\ H k+i^ n) + \\tpr\\- 



Using local Poincare inequalities and the stability (4.3) of ip r we get 



\<Pr 



<||W^||<^ m 



i\\L 2 (Fi) 



/N P . \ 2 

= hi \ J2 ll Vtt '" - Vi sz u ■ nWl^Fi)) ■ 

Using the stability of the projection onto piece wise constants, element wise trace 
inequalities and finally approximation, we conclude 



Vw • n — Vig Z it • n 



\h {Fi) < || Vu • n - Vi k z u ■ n\\ 2 L 2 iFi) 



< 2CUh- 1 \\V(u-t k sz u)\\ 2 L2(Pi) + h ]T \\D 2 (u-i^u)\\l 2{K) ) < ft 2fc -VlW(p ( ) 



KePi 
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where D 2 u is the standard multi-index notation for all the second derivatives of u. 
We conclude that 



\tPrW ^ ($2 \\Vu ■ n - Vig Z u ' Hli 2 ^))^ ~ h k+1 \u\ H k+ 



The estimate on the gradient is immediate by 



|V(u - n a u)\\ < ||V(u - &u)|| + l|V(4« - t»«)II 

< ||V(u - «sz u )ll + C/ /l_1 ||isz M - ^"11 < ^Mfffe 



i(f2)- 



□ 



6.2.1. Convection dominated regime: the Continuous Interior Penalty 
mehod. In this section we will sketch how the above results extend to symmetric 
stabilization methods assuming that c a > 0. To reduce technicalities we also assume 
that ft £ M 2 . We give a full proof only in the case of piecewise affine finite elements. 
Recall that the CIP method is obtained by adding a penalty term on the jump of the 



gradient over element faces to the finite element formulation (6.2). The formulation 
then writes: find Uh G V k such that 



A h (u h ,v h ) + Jh(u h ,v h ) = (f,v h )n, Vv h e V£, 



(6.12) 



where 



Jh{uh,v h ) 



KeT h F€dK\dQ 



h 2 F \(3 ■ n F \\S7u h ■ n F ][Vv h ■ n F ] ds, 



with [x] denoting the jump of the quantity x over the face F and np the normal to 
F, the orientation is arbitrary but fixed in both cases. 

The analysis once again depends on the construction of a special interpolant 
kcipu G Vh- This time ttqipu must satisfy both the optimal approximation error 



design condition: 



estimates of (6.9), the property (6.10) on the normal gradient, and the additional 



(u - ttcipu, /3 ■ Vv h ) < \\h 



: {u-Tr C ipu)\\7 CIP J h (v h ,v h )2, Wv h G V h . 



(6.13) 



Once such an interpolant has been proven to exist, the technique of [3J, combined with 
the analysis above, may be used to prove quasi-optimal L 2 -convergence for c a > 0. 
Using a similarly designed interpolation operator, an inf-sup condition can be used to 
prove stability and error estimates in the norm ||| • \\\h,s following [6j (5] . Here we will 
first the error estimate in the L 2 -norm, assuming the existence of Ttcipu and then 
show how to construct the interpolant in the special case k = 1. 

Proposition 6.6 
exists. Let u G H +1 (f2) be the solution to (6.1), with c 
to (16.121). Then 



Assume that ttqipu G V k satisfying, (6.9), (6.10) and (6.13) 



> 0, and Uh be the solution 



\\u-u h \\ < ( C(T )-\aihi +|#(1 + Pe-^))h k+ i\u\ Hk+HQ) . 
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Proof. Let e/,. := Uh — ncipu. There holds with c„ > 0, 

Co- He/,,11 2 + IIKIH^o + J h (e h ,e h ) < A h (e h ,e h ) + Jh{e h , e h ). 
By the consistency of the method we have 

Calle^H 2 + \\\e h \\\ 2 h o + Jh(eh,eh) < A h {u - ncipu,e h ) - J h {^cipu,e h ). 



Finally by the continuity (6.11), that holds thanks to property (6.10), we have 



A h (u - itcipu, e h ) - J h (ir C ipu, e h ) 

= (cr(u - ir C ipu),e h ) + (u - wcipu, (3 ■ Ve h ) - I f3 ■ n(u - n C ipu)e h ds 

Jan 

+ ea h (u - ncipu, e h ) + J h (Tr C ipu, e h ) 
< {(a*h? + C\f3\^j~f p )\\h-^ (u - tt C ipu)\\ + \\u - TTcipu\\i, h ,p 
+ e^\\u - 7r C /pu||* + Jhi^cipu^cipu) 1 ) 

x (a^\\e h \\ 2 + \\e h \\l h ^ + J h (e h ,e h ))i (6.14) 

and we end the proof by applying approximation estimates. □ 

We will now prove the existence of the interpolant ncipu in the case of piecewise 
affinc continuous finite element approximation. 



Lemma 6.7. The function ire ip u S V^, satisfying (6.9), (6.10) and (6.13), is 
well defined and satisfies the approximation estimate 

\\u - Tr C ipu\\ + h\\V(u - -Kcipu)\\ < h 2 \u\ H 2 (n y 



Proof. We write ttqipu :— tt^u + fcip where tt^u denotes the L 2 -projection on 
and (pcip € is a function defined on patches Pi that satisfy the inequalities 
(4.2) and (4.3), but also has the property 



ipcip dx = 0, i = 1, . . . , N P . 



Clearly for this to hold we must modify the definition of the patches on the faces Fi 
to include interior nodes in the domain. For simplicity we assume that any element 
containing a node that connects to two nodes in the boundary segment Fi (through 
edges that may be associated to other elements) is included in the patch Pi (see Figure 
6.1). Define two functions wj and Wp on Pi (also illustrated in Figure [6~T|) such that 



Wl 



1 in all nodes x € P, 

o 

in all nodes x € O \ P, 



Wp 



1 in all nodes x € Fi 

_ o 

in all nodes x € ft \ Fi 



We must now show that there exists a function (fi — awi + bwp satisfying the two 
constraints 



ipi dx = 0, Vtp l ■ n = rv 



(6.15) 



The construction of ncipu is obtained by choosing = V« ■ n — V7T/jW • n in the 
system (6.15) above and then defining fcip\pt '■— fi- 
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To study ifi, first map the patch Pj to the reference patch Pj. Consider the linear 
system for v := (a, b) T £ R 2 of the form: 



Av 



J p wi dx 



dx 




a 


n ds 




b 



f F V(u - KhU) ■ n ds 



=:/• 



We must prove that the matrix A is invertible, but this is immediate noting that the 
two coefficients in the first line of the matrix both are strictly positive, whereas in the 
second line the coefficient in the first column is negative by construction and that in 
the right column is positive. The stability estimate (4.3 1 now follows from a scaling 
argument back to the physical patch Pj. Indeed since the matrix A is invertible we 
have 



< 



sup 



w Av 



sup 



\W\ „. ' 2 if 

By norm equivalence we have 

\\0iWp, £ l|V0i||p. 
After scaling back to the physical element we get 



l/l- 



< 



< 



l/l- 



< 



\Vtp, 



i\\Pi 



< 



|/| < \\h*V(u-Tr h u)-n\\ Fi , 



(6.16) 



which proves (4.3 1. 

The approximation error estimates are proven in the same way as in Lemma |6.5 
Indeed by a similar decomposition of the error we have for this case 



■u. ■ 



■kcipu\\ < \\ u - nh,u\\ + \\-KhU - ircipu\\ < h 2 \u\ H 2^ + \\(ficip\\ 



and for <ficip we may conclude using the proof of Lemma 6.5 using (6.16). 



It remains to prove the continuity (6.13). This follows from 
(u - ncipu, p ■ Wv h ) = (u- ir h u, @ ■ Wv h ) + ^(<ft, P ■ V« fc ) 



i=i 



N P 



(u - n h u, p ■ Wv h - IcipP ■ Vv ft ) + y^(y>i, (P ■ Vu h - 7T 0: p./3 • VVfc)). 



1=1 



Here Iqip denotes a particular quasi interpolation operator defined using averages of 
P ■ Vw/j in each node (see [3]) and iro,P t denotes the projection on piecewise constant 
functions on Pj. Using norm equivalence on discrete spaces and mapping from the 
reference patch, we observe that 



and 



N P 

5>*i 

i=l 



! -(P ■ Vv h - IcipP ■ Vv h )f < -ycip J h( v h,v h ) 



■5(/3- Vv h -TC 0tPi P- Vv h )f P . 
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< 



lC!P J h(Vh,Vh)- 



Fig. 6.1. Example of a boundary patch Pi, with the functions wj (left) and uip (right). The 
functions take the value 1 in filled nodes and zero in the other nodes. 



N 


Nitsche H 1 


strong H 1 


Nitsche L 2 


strong L 2 


10 


7.0E-1 (— ) 


6.7E-1 (— ) 


2.4E-2 (— ) 


2.0E-2 (— ) 


20 


3.5E-1 (1.0) 


3.5E-1 (0.94) 


5.5E-3 (2.1) 


5.5E-3 (1.9) 


40 


1.7E-1 (1.0) 


1.7E-1 (1.0) 


1.3E-3 (2.1) 


1.3E-3 (2.1) 


80 


8.2E-2 (1.1) 


8.2E-2 (1.1) 


3.3E-4 (2.0) 


3.1E-4 (2.1) 



Table 7.1 

Comparison of errors between the non-symmetric Nitsche method and standard strongly imposed 
boundary conditions, using piece wise affine approximation on unstructured meshes. 



The first claim was proved in [3] and the second holds since /3 • Vw^ is constant on 
each element. □ 

Remark 1. For high order element the construction of the interpolant ttqipu is 
much more technical and beyond the scope of the present work. Indeed it is no longer 
sufficient to prove orthogonality of ipi against a constant on Pi, but it must be shown 
to be orthogonal to the continuous finite element space of order k — 1 , on Pi. On the 
other hand the patches Pi can be chosen freely, provided diam(Pi) = 0(h). 

7. Numerical examples. We study two different numerical examples, both 
have been computed using the package FreeFem++ 12). First we consider a simple 
problem with smooth exact solution, then we consider a convection-diffusion prob- 
lem and show the stabilizing effect of the Nitsche type weak boundary condition for 
convection dominated flow. 



7.1. Problem with smooth solution. We consider equation (2.1) in the unit 
square, with / = 5ir 2 sin(7r:c) sin(27ry) and g = 0. The mesh is unstructured with 
N = 10, 20, 40, 80 elements per side. The exact solution is then given by u = 
sin(7ra;) sin(27ry). We give the convergence in both the L 2 -norm and the i/^-norm 
for piece wise affine approximation in Table [7Tj The case of quadratic approximation 



is considered in Table 7.2 The order p in 0(h p ) is given in parenthesis next to the 
error. We have not managed to construct an example exhibiting the suboptimal 
convergence order of the Nitsche method. Some cases with non-homogeneous bound- 
ary conditions, not reported here, were computed both with affine and quadratic 
elements. They all had optimal convergence on the finer meshes. The theoretical 
results do not extend to the symmetric version of Nitsche's method and stability is 
unlikely to hold on general meshes. Applying the symmetric method to the proposed 
numerical example yields a solution with clear boundary oscillations on the coarse 
meshes see Figure |7.1| On finer meshes these oscillations vanish and the performance 
is similar to that of the non-symmetric method. Note that although the convergence 
of the Nitsche method is optimal in this case, the error constant of the non-symmetric 
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N 


Nitsche H 1 


strong H 1 


Nitsche L 2 


strong L 2 


10 


5.3E-2 (— ) 


5.1E-2 (— ) 


1.7E-3 (— ) 


6.5E-4 (— ) 


20 


1.4E-2 (1.9) 


1.4E-2 (1.9) 


2.2E-4 (2.9) 


9.6E-5 (2.8) 


40 


3.5E-3 (2.0) 


3.5E-3 (2.0) 


2.1E-5 (3.4) 


1.1E-5 (3.1) 


80 


8.6E-4 (2.0) 


8.6E-4 (2.0) 


2.5E-6 (3.1) 


1.4E-6 (3.0) 



Table 7.2 

Comparison of errors between the non- symmetric Nitsche method and standard strongly imposed 
boundary conditions, using piece wise quadratic approximation on unstructured meshes. 



error norm 


7 = 


7 = 10 


7 = 20 


7 = 40 


7 = 80 


\\u - U h \\ L 2 


3.3E-4 


2.9E-4 


3.0E-4 


3.0E-4 


3.0E-4 


\\u - U h \\ H l 


8.2E-2 


8.2E-2 


8.2E-2 


8.2E-2 


8.2E-2 



Table 7.3 

Study of the dependence of the accuracy on the penalty parameter, piece wise affine approxi- 
mation, unstructured mesh, N = 80 



method in the L 2 -norm is a factor two larger than that of the strongly imposed bound- 
ary conditions for piece wise quadratic approximation. The same computations were 
made on structured meshes (not reported here) and this effect was slightly larger in 
this case, with a factor two in the affine case and four in the quadratic case. The 
errors in the H ^norm on the other hand are of comparable size for the two methods. 
This motivates a study of how the error depends on the penalty parameter 7 in 



(3.3). We therefore run a series of computations with 7 = 0,10,20,40,80. In Table 



7.3 we report the results for piece wise affine approximation and in Table 7.4 the 
results for piece wise quadratic approximation. We note that there is a visible, but 
negligible, effect on the error measured in the L 2 -norm, but no effect on the error in 
the _ff 1 -norm. 

7.2. Problem with outflow layer. For this case we only compare the solutions 



qualitatively. We consider the problem with a convection term (6.1 1. To create an 
outflow layer we have chosen / := 1, :— (0.5, 1), a := in ft. We discretized ft with 
a structured mesh having 80 piece wise affine elements on each side. The contour- 
plots for e — 0.1, 0.001, 0.00001 are reported in Figure [73] for Nitsche's method and in 
Figure |7.3| for the strongly imposed boundary conditions. Note that no stabilization 
has been added in either case. This computation illustrates the strong stabilizing 
effect of the weakly imposed boundary condition. A theoretical explanation of this 
phenomenon was given in [18 . Finally we consider the effect of adding stabiliza- 
tion to the computation. In this case we take N = 80 with piece wise quadratic 
approximation. We report the results of a computation without stabilization, with 



the SD-method (<y SD = 0.2) and with the CIP-method (70/p = 0.005) in Figure [7^4 
Note that the stabilized methods clean up the remaining spurious oscillations in both 
cases. 

Acknowledgements This note would not have been written without Professor Tom 
Hughes who told me that the non-symmetric Nitsche's method appeared to be stable 
without penalty in large-eddy simulations and pointed me to the reference [13] . I 
would also like to thank Professor Rolf Stenberg for interesting discussions on the 
subject of Nitsche's method. 



19 



error norm 


7 = 


7 = 10 


7 = 20 


7 = 40 


7 = 80 


\\u - Uh\\L 2 


2.1E-5 


1.3E-5 


1.2E-5 


1.2E-5 


1.2E-5 


\\u - U h \\ H l 


3.5E-3 


3.5E-3 


3.5E-3 


3.5E-3 


3.5E-3 



Table 7.4 

Study of the dependence of the accuracy on the penalty parameter, piece wise quadratic approx- 
imation, unstructured mesh, N = 40 




Fig. 7.1. Comparison of the contourplots of the unstabilized non-symmetric method (left) and 
symmetric (right) method, piece wise affine approximation, N=10. 




Fig. 7.2. Convection- diffusion equation discretized using the non-symmetric Nitsche boundary 
condition, N = 80, piece wise affine approximation, from left to right: e = 0.1, e = 0.001, e = 
0.00001. 
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Fig. 7.4. Convection- diffusion equation discretized using the non- symmetric Nitsche bound- 
ary condition, N = 80, e = 0.00001, piece wise quadratic approximation, from left to right: no 
stabilization, SD- stabilization (-fsD = 0.5j, CIP- stabilization (-fciP = 0.005j. 
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